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>^ Abstract 

c3 : 

This is a chapter of the multi-author book "Understanding Quantum Phase Transitions," edited 
by Lincoln Carr and published by Taylor & Francis. In this chapter, we give a general introduction 
to the worm algorithm and present important results highlighting the power of the approach. 

O 

<D • PACS numbers: 

B' 

> 

a 

a 

g 

o 

m 
> 
m 

On 

m 



ON 

o 



X 



i 



Theoretical studies often involve mappings of the original system onto an equivalent (with 
regards to the final answer for some property) description in terms of abstract mathemati- 
cal/graphical objects. Path integrals, high-temperature expansions, and Feynman diagrams 
are the well-known examples considered in this chapter. Under mapping, one has to deal 
with the infinite-dimensional configuration space having complex topology and non-local 
constraints which severely reduce efficiency of Monte Carlo simulations based on standard 
local updates. This sometimes leads to ergodicity problems in large system when the entire 
configuration space can not be sampled in a reasonable computation time. A somewhat re- 
lated difficulty facing conventional Monte Carlo schemes is the computation of off-diagonal 
correlation functions since they have no direct relation to the configuration space of the 
partition function. In what follows we consider path integrals for lattice and continuous 
systems, high-temperature expansions, and Feynman diagrams, explain the general idea 
of how Worm Algorithms deal with the topological constraints by going to the enlarged 
configuration space, and present illustrative results for several physics problems. 



PATH-INTEGRALS IN DISCRETE AND CONTINUOUS SPACE 

For clarity, we start by introducing the Hamiltonian describing lattice bosons (a straight- 
forward generalization to fermions will be mentioned later) making hopping transitions be- 
tween the nearest-neighbor sites < ij >, interacting by the pairwise potential Uik, and 
subject to the external potential /if 

H = + H' = ±J2 U ^ n ^ - E W ~ t E h h ■ (!) 

i,k i <ij> 

Here bi is the bosonic annihilation operator and n« = b\bi. In the Fock basis of site occupation 
numbers, \ct) = the first two terms, representing H^°\ are diagonal, while the last 

term, representing H', is not. We write the statistical operator as 

c -/jh = c -/jh(o) exp r _ r fi dTH , {T) ) (2) 







where H'{t) = e TH(0) H'e tH(0 \ and the exponential is understood as the time-ordered ex- 
pansion 

1 - / dr H'(t) + / dn dr 2 H'lrjH'fo) + ■■■ ■ (3) 

JO JO Jt! 



Lattice path-integrals immediately follow from the graphical representation of the expan- 
sion ([3]). Consider the partition function Z given by the trace of e~^ H and write explicitly 
the m-th term (before integration over time) as 

(-I)™ d-r e -(^X (^ o a J e -(n-. 2 )< (tf^J e -^> , (4 ) 

with a m = «o to reflect the trace condition (periodic boundary condition in imaginary time), 
and d m T = dr\ . . . dr m . Since hopping terms in H' change the state by shifting only one 
particle to a nearest-neighbor site, the sequence of matrix elements is completely determined 
by specifying the "evolution" or "imaginary time trajectory" of occupation numbers {n^r)}. 
In the left panel of Fig. [T]we show the trajectory describing one of the 4-th order terms which 
contributes t 4 d 4 r 1 ■ 2 ■ \pl ■ a/2 exp{— drH^(r)} to Z, where for brevity we use H^(r) 
for the energy of the state |{^(t)}) and give explicit expressions for the hopping matrix 
elements (n* — 1, rij + 1| — tb^b^rii, rij) = —ty/n^rij + 1). Thus, the partition function can 
be written as a sum over all possible paths {rij(r)} such that rij(/3) = rij(O) 

Z= £ W[{ ni (T)}}, (5) 

K(t)} 

with strict rules relating the trajectory shape to its contribution to Z. The trajectory weight 
is sign-positive if t is positive or the lattice is bi-partite. 

The path-integral language at this point is nothing but a convenient way of visualizing 
each term in the perturbative expansion ()3]). Due to the particle number conservation, 
the many body trajectory can be decomposed into the set of closed (in the time direction) 
single-particle trajectories, or worldlines. Worldlines can "wind" around the /3-circle several 
times before closing on themselves. In a system with periodic boundary conditions in space, 
the trajectories can also wind in the space direction. Worldlines with non-zero j3- and 
space- winding numbers are said to form exchange cycles and winding numbers, respectively; 
they are directly responsible for superfluid properties of the system [l| and are the origin of 
non-local topological constraints mentioned in the introductory paragraph. 

The Green's function of the system, 

G{t m - 77, i M - ii) = r r (6j M (r M )& ir (T/)) , (6) 

has a similar path-integral representation, see right panel in Fig. (TJ with one notable differ- 
ence: due to operators at point {im-i'Tm) and b at point (ij, tj) there is one more particle 
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FIG. 1: Lattice path-integral representations for the partition function and the Green's function. 
Line thickness is proportional to rij 

present in the system on the time interval (17, tm), i- e. when the {^(t)} evolution is 
decomposed into worldlines there will be one open worldline originating at point tm and 
terminating at point 77. These special points will be labeled as X ("Ira") and M. ("Masha") 
throughout the text. We will use short-hand notations, Z-path, and G-path to distinguish 
configuration spaces of Z and G. 

We skip here the standard derivation of Feynman's path-integral representation for an 
interacting iV-particle Hamiltonian in continuous space 2] 

1 JV 1 JV N 
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leading to 



VH T exp 



mR 2 



+ U(R) - fi(R) 



dr 



(8) 



where R T = {77 (r)} represents positions of all particles at time r, U(R) and —fi(R) stand 
for internal and external potential energy respectively, and Rp = i?o to satisfy the trace 
condition. In practice, the imaginary time axis is sliced into sufficiently large number of 
intervals and the trajectory is defined by specifying particle positions at a discreet set of 
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FIG. 2: Loop representations for the partition function and the correlation function for the Ising 
model. Line thickness is proportional to N(,. 

time points. Apart from time slicing, there is no fundamental difference between the lattice 
and continuous path-integrals. The Z-path consists of closed worldlines (the decomposition 
of the many-body path into individual worldlines is unique in continuous space), and the 
G-path contains one open worldline originating at (tm,tm) and terminating at (77, tj). 

The only difference between the fermioninc and bosonic systems is in the sign rule: for 
fermions the trajectory weight W involves an additional factor (— l) p , where p is the parity 
of the permutation between the particle coordinates in Rp relative to Rq. 



LOOP REPRESENTATIONS FOR CLASSICAL HIGH-TEMPERATURE EXPAN- 
SIONS 

Classical statistical models can be also mapped to the configuration space of closed loops. 
Important examples include Ising, XY, and \ip\ A models, as well as their multi-component 
generalizations with or without the gauge coupling, e.g. the CP 1 model. We illustrate 
the general idea of the mapping by considering the simplest case of the Ising model when 
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—H/T = K^2 <i j > a i crj for N Ising spin variables cij = ±1. One starts with the partition 
function and expands bond Gibbs factors into Taylor series (to simplify notations we will 
use subscript b for lattice bonds) 

z = E n E n eCW 1 - o»> 

{<Tj=±l} b=<ij> {o- l= ±l} b=<ij> N b =0 

By changing summation over {iV;,} and {<7j = ±1} places we obtain 

loops tsN loops 

Z = 2 N T, II JfT = 2 * £ W H N b}} ■ (10) 

{N b } b=<ij> b ' {N b } 

The "loops" label on the sum represents the constraint that the sum of all bond numbers in- 
cident on every lattice site, Li = J2b=<ij> Nb, nas ^° ^ e even i otherwise Y^ ai ±i aLi ls zero - 111 
the graphical representation where iVj, is substituted with Ni, lines drawn on the correspond- 
ing bond, this is equivalent to demanding that the allowed configuration of lines is that of 
closed un-oriented loops, since loops always contribute an even number to Lj = J2b=<ij> ^b, 
see the left panel in Fig. HJ [Using identity e Kai(j3 = cosh(i^) J2N b =o i [tanl^iQcrjcrj]^ 6 we 
get a more compact formulation since only one line can be drawn on the bond.] In close 
analogy with the Green's function, the configuration space of the spin-spin correlation func- 
tion Gjm = (o'jO'm) is that of closed loops with one open line originating at site I and 
terminating at A4, see the right panel in Fig. [2j The difference between the XY and Ising 
models is that loops are oriented in the XY-case, i.e. iVj G (— oo, oo). 



WORM ALGORITHM: THE CONCEPT AND REALIZATIONS 

The Worm Algorithm (WA) strategy for updating loop configuration spaces involves two 
major ideas: 

1. The configuration space is enlarged to include one open line, as if someone started 
drawing a new loop but is not finished yet. In all examples mentioned above this is not 
merely an algorithmic trick but also an important tool to have direct access to off-diagonal 
correlation functions, such as the Green's function j^, 0] (with two worms one may calculate 
multi-particle off-diagonal correlations as well jsl, 0].) From time to time the two ends of the 
open line come close and get connected thus making a loop and transforming G-path to Z- 
path. In Monte Carlo, drawing and erasing processes are balanced and are complementary 
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to each other. 

2. All updates on G-paths are performed exclusively through the end-points of the open 
line, no global updates or local updates transforming one Z-path to another Z-path are 
necessary. 

More generally, the Worm Algorithm idea is to consider an enlarged configuration space 
which includes structures violating constraints present in the Z-sector of the space. Of- 
ten, this can be achieved automatically by considering the relevant correlation functions. 
Green's function and, for paired states, higher-order off-diagonal correlators are the ap- 
propriate choice for the path-integral space. However, one should feel free to introduce 
"unphysical" configurations which bear no meaning at all and are used solely for employing 
local updates to produce non-trivial global changes of physical configurations. An example 
with the momentum conservation law in Feynman diagrams discussed below illustrates the 
point. 

WA is a local Metropolis scheme, but, remarkably, its efficiency is similar to (or better 
than) the best cluster methods at the critical point, i.e. it does not suffer from the critical 
slowing down problem. It has no problem producing loops winding around the system, al- 
lows efficient simulations of off-diagonal correlations, grand canonical ensembles, disordered 
systems, etc. Below we provide more specific details of how it works. 

Discrete Configuration Space: Classical High- Temperature Expansions 

We start with the simplest case of the classical Ising model. The entire algorithm consists 
of just one update: 

If X = Ai, select at random a new lattice site j and assign X = Ai = j; otherwise skip 
this step. [In other words, put your pensil/eraser anywhere.] Select at random the direction 
(bond) to shift Masha to a n.n. site, let it be k, and propose to change the bond number from 
N b to N' h = mod 2 (A" 6 + 1). Accept the move with probability R = max 
This is a complete description of the algorithm! 

Every configuration contributes unity to the statistics of G^m- For the Ising model 
Z = Gi=m- Generalizations to other classical statistical models are straightforward 
Ergodicity is guaranteed because a finite number of steps is required to erase any initial 
trajectory and then to draw a new one, line after line. The efficiency is ultimately linked to 
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FIG. 3: The upper and lower panels illustrate transformations performed by the Open/Close and 
Insert /Remove pairs of updates respectively on the lattice path-integral configurations. 

the fact that WA works directly with the correlation function of the order parameter field. 
Continuous Time: Quantum Lattice Systems 

We now illustrate the WA updating strategy for the system of lattice bosons described 
by Eqs. ([I]), fl5]) with the configuration space shown in Fig. HJ The set of updates presented 
below forms an ergodic set. It is sufficient to describe updates performed with the Masha- 
end of the open worldline; updates involving the Ira-end follow immediately from the time 
reversal symmetry. 

Open/Close. — This pair of updates takes us back and forth between the Z- and G-paths 
by selecting an existing worldline and erasing a small part of it (open) or drawing a small 
piece of worldline between the end-points of the open line to complete the loop (close). 
These updates are illustrated in the upper panel of Fig. [3] In the (open) update the path 
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FIG. 4: The upper and lower panels illustrate transformations performed by the Move and 
Jump /anti-Jump updates respectively on the lattice path-integral configurations. 

interval characterized by time-independent occupation numbers on a given site is selected 
at random from the list of such intervals, and the imaginary times tj < tm are seeded 
from the normalized probability distribution p(rj, tm) with n, tm G (r in , rg n ). Formally, the 
distribution function is arbitrary, and this freedom should be used to optimize the acceptance 
ratio R oc W n ew/\W idP(Ti, tm)\- in the {close) update one checks whether Ira and Masha 
are connected by a single path interval with time-independent occupation numbers on a 
given site, and if true, proposes to eliminate them by "connecting the dots" . We skip here 
further technical details (which are minimal), as well as explicit expressions for acceptance 
ratios, which can be found elsewhere ^]. 

Insert/Remove . — This pair of updates also switches back and forth between the Z- and 
G-paths by drawing a small piece of a new worldline {insert) or erasing a small piece of 
worldline between the end points {remove). The two updates are illustrated in Fig. [3] and 
in all respects are similar to the Open/Close pair except that the time ordering of 77 and 
Tm is reversed, see the lower panel in Fig. |3j 
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FIG. 5: The upper and lower panels illustrate transformations performed by the Open/Close and 
Insert/Remove pairs of updates respectively on the continuous-space path-integral configurations. 

Move. — Once in the G-path space the algorithm is proposing to move Masha along the 
time axis, tm — > t' m within the bounds determined by the change of the occupation numbers 
on a given site. This motion is equivalent to drawing/erasing the worldline, see the upper 
panel in Fig. HI 

Jump /anti- Jump. — This is the only complementary pair of updates which is different 
from continuous transformations of lines since it involves the motion of end-points in space. 
Without changing the time position of Masha, we place it on the neighboring site and 
connect worldlines of the two sites involved in the update in such a way that the rest of the 
path remains intact. This requires adding/removing a kink immediately before or after r^/. 
We illustrate the Jump /anti- Jump updates in the lower panel of Fig. HJ Note the difference 
between the two cases. Whenthe kink is inserted to the left of Masha, the transformation 
can still be interpreted as proceeding with drawing the same worldline. When the kink is 
inserted to the right of Masha, we reconnect existing worldlines and ultimately effectively 
sample all allowed topologies of the path. 
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FIG. 6: The upper and lower panels illustrate transformations performed by the Draw/Erase and 
Swap updates respectively on the continuous-space path-integral configurations. 

This concludes the description of the algorithm. Such properties as density, energy, 
density- density correlations, etc. are computed using standard rules when the configuration 
is in the Z-path sector. Every G-configuration makes a direct contribution to the statistics 
of G(t m - T T ,i M - ii)- 

Bosons in Continuous Space 

Worm Algorithms for the continuous and lattice systems are essentially identical at the 
conceptual level. Most differences are technical and originate from having discreet instead 
of continuous time and continuous instead of discrete space. Specific protocols of how one 
proposes new variables and performs measurements can be found in Ref. [7J. The only dif- 
ference worth mentioning is that in continuous space the decomposition of the path into 
individual worldlines is unique (which actually simplifies things). Here we simply illustrate 
the updates. Since in all cases we know how to compute the path contribution to the statis- 
tics of Z or G, graphical representations can always be converted to precise mathematical 
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FIG. 7: Effective mass for particle (circles) and hole (squares) excitations in the 2D Bose Hubbard 
model at unity filling as a function of the hopping-to-interaction ratio J/U. The exact results at 
J/U = are m+ = 0.25/ J and m_ = 0.5/ J. The dashed lines show the lowest-order in J/U 
corrections to the effective masses. Close to the critical point the two curves overlap, revealing the 
emergence of the particle-hole symmetry implied by criticality-induced emergent Lorentz invari- 
ance. The sound velocity of the relativistic spectrum in the Lorentz-invariant regime is found to 



be c/J = 4.8 db 0.2. (Reproduced from Ref. 



expressions for the acceptance ratios which account for the ratio of the configuration weights 
and probabilities/probability densities of applying a particular type of update. 

Figures |5] and [6] show an ergodic set of updates which would allow one to efficiently sim- 
ulate continuous space systems. The Open/Close and Insert/Remove pairs of updates are 
reminiscent of those in lattice systems. The Draw/Erase pair naturally combines in one up- 
date both space and time motion of the end-point and is essentially a literal implementation 
of the draw-and-erase algorithm. The Swap update is equivalent to the version of the Jump 
update which involves reconnection of the particle worldlines. In Swap all modifications of 
the path are restricted to occur between the two dashed lines, see Fig. [6] 

The rules for collecting statistics to the diagonal and off-diagonal properties are similar 
to lattice models. 



Momentum Conservation in Feynman Diagrams 



Diagrammatic Monte Carlo (see, e.g., Ref. [8| and references therein) is a technique of 
sampling entities expressed in terms of Feynman's diagrammatic series by a Markov process. 
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The configurational space of the process consists of Feynman's diagrams with fixed internal 
variables. — It is the sampling process itself that accounts for the integration over these vari- 
ables, on equal footing with summation over the order and topology of the diagrams. In the 
process of producing the Markov chain of diagrams, generating the (n+l)-st diagram is per- 
formed by applying one of a few elementary updates to the n-th diagram. The elementary up- 
dates are supposed to change the structure of a diagram (by adding/removing/reconnecting 
a small number of propagators) and the values of (a small number of) internal variables. 
In the space-time representation, updating internal variables causes no problem since these 
are nothing but the space-time points corresponding to the ends of propagators, their val- 
ues being naturally generated with appearance of new propagators, and abandoned when 
the corresponding propagators are removed. In momentum and/or frequency representa- 
tion, the situation is quite different. The momentum (for briefness, from now on we speak 
of momentum only) of a given propagator is not independent of the momenta of other 
propagators in view of the momentum conservation constraint taking place at each vertex. 
Adding/removing a propagator to/from a diagram would thus involve a change of momenta 
of other propagators, rendering the updating routine complicated and less efficient. 

A simple and efficient way out is provided by the "worm" idea of working in an extended 
configurational space. The additional class of diagrams that we need consists of diagrams 
featuring worms, by which here we mean vertices with non-concerving momenta. Clearly, 
the minimal non-trivial number of worms is two. Normally, working with no more than two 
worms proves sufficient. Note that if the algebraic sum of all the momenta at one of the two 
worms is 5, then its counterpart at the other worm is —5, so that S is the only continuous 
parameter associated with the pair of worms. 

A crucial observation now is that if all the structural updates of diagrams are per- 
formed in the subspace of diagrams with worms and the ends of the propagator (s) to be 
added/removed/reconnected are linked to the two worms, then the worms will readily "ab- 
sorb" the residual momentum, k res , associated with the update, the only consequence for 
the worms being 5 — » 5 + k vcs . Details on implementing this idea can be found in Ref. js]. 
The overall updating scenario is as follows. Switching between physical and worm sectors of 
the configurational space is achieved by a pair of complementary updates creating/deleting 
a pair of worms at the ends of a propagator, with simultaneously changing the momentum 
of this propagator. The rest of the updates is performed in the worm sector. The efficiency 
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of the scheme is achieved by introducing an update that translates worms along one of the 
propagators attached to the worm vertex. In the updated diagram, the conservation of mo- 
mentum at the worm's original position is ensured by appropriately changing the momentum 
of the propagator along which the translation is being performed. Translating worms along 
propagators allows one to efficiently sample all their positions and thus apply the updates 
changing the structure of the diagram — associated with the worms, as discussed above — in 
a generic way. 



ILLUSTRATIVE APPLICATIONS 
Optical-Lattice Bosonic Systems 

Bosons in optical lattices, being accurately described by the Hubbard Hamiltonian (CQ) 
9|, are perfectly suitable for simulations by worm algorithm. With a standard desktop 
computer, one can simulate equilibrium properties of 3D systems with 200 3 lattice sites. In 
Refs. 
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ll| it has been demonstrated that this approach allows one to obtain precision 
results for equations of state and produce an accurate phase diagram of the system. It is also 
possible to trace the evolution of the particle/hole spectrum of elementary excitations with 
decreasing interaction strength, from the strong coupling limit down to the critical point 
of the Mott-insulator-to-superffuid quantum phase transition. Especially interesting is the 
vicinity of the quantum critical point, where the emergent Lorentz invariance brings about 
particle-hole symmetry, see Fig. [7J Along with the pure single- component bosonic Hubbard 
model, one can simulate mult i- component and disordered systems, see Figs. [3191 

The lattice bosonic systems are in the focus of Optical Lattice Emulator project supported 
by DARPA and aimed at the development, within the next few years, of experimental tools 
of accurately mapping phase diagrams of lattice systems by emulating them with ultracold 
atoms in optical lattices. The numerically exact solutions for real experimental systems 
will be used for validating the emulators. The first successful validation of the emulator 



of the Bose Hubbard model was reported in Ref. [13|. At the heart of the protocol is the 
direct comparison of the experimental time-of-ffight images with the theoretical ones. The 
latter are produced by time-evolving the single-particle density matrix obtained in a direct 
simulation of a given number of atoms in a trap, see Fig. [101 
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FIG. 8: Groundstate phase diagram of 2D two-component bosonic Hubbard model at half-integer 
filling factor for each component. The on-site interactions within each component are infinitely 
strong (hard-core limit), while the inter-component interaction, U, is finite. The hopping elements 
of the two components are t a and tb, the parameter z = 4 is the coordination number. The revealed 
phases are as follows, (i) checkerboard solid in both components, a.k.a. z-Neel phase (2CB), (ii) 
checkerboard solid in one component and superfluid in its counterpart (CB+SF), (hi) superfluid 
in both components (2SF), (iv) super-counter-fluid, a.k.a. Xy-ferromagnet (SCF). The observed 
transition lines are: 2CB-SCF (first-order), SCF-2SF (second-order), 2CB-2SF (first-order), 2CB- 
CB+SF (second-order), and CB+SF-2SF (first-order). Lines are used to guide an eye. (Reproduced 
from Ref. [6|.) 

Supersolidity of Helium-4 

A supersolid is a quantum solid that can support a dissipationless flow of its own atoms. 
Here the term 'solid' is understood in the most general context of any (regular or amorphous, 
continuous-space or lattice) state with broken translation symmetry. 
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FIG. 9: Groundstate phase diagram of the disordered 3D Bose-Hubbard model (at unity filling) in 
the vicinity of the point of the superfiuid(SF)-to-Mott-insulator(MI) quantum phase transition. 
The E g / 2 (U) curve marks the Bose-glass(BG)-Mott-insulator transition boundary according to the 
conjecture that the transition occurs when the bound of disorder reaches the half-gap, E g / 2 , of the 
pure Mott insulator. Error bars are shown, but are smaller than point sizes. (Reproduced from 
Ref. Irj.) 



The modern age of supersolidity of bosonic crystals in continuous space began with 
;he discovery by Kim and Chan of non-classical rotational inertia (NCRI) in solid 4 He 
llo] ]; for reviews of further developments in the field, including direct observation of 
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and references therein. In the 



a superflow, as well as preceding work, see Refs. 
combined experimental and theoretical effort aimed at understanding the microscopic picture 
behind the effect, the first-principles simulations of regular and disordered solid 4 He play 
a very important part. Here we present some numeric results shedding a direct light on 
microscopic mechanisms of supersolidity in 4 He. 

To introduce a general theoretical background for interpreting numeric results, we start 
with a number of rigorous statements. From the field-theoretical perspective, a super- 
fluid/supersolid groundstate — as opposed to an insulating one — is almost trivial, since the 
topological constant of motion responsible for the superfluidity is naturally introduced in 
terms of the phase of the classical matter field, so that phenomenon of superfluidity in a 
quantum field is simply inherited from the classical counterpart of the latter. An insulating 
groundstate is possible only in a quantum field. It is an essentially non-perterbative and 
thus non-trivial phenomenon; the fact of its existence in bosonic systems has been recently 
denied by P.W. Anderson [? ]. 
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FIG. 10: Comparison of experimental and simulated time-of-flight distributions: Shown the inte- 
grated column density n±(x, y) represented by the optical density as obtained from the experiment 
and the QMC simulations for different temperatures and two lattice depths. (Reproduced from 
Ref. |13( | : see this reference for more detail.) 

A proof of existence of insulating groundstates in bosonic crystals immediately follows 



from the theorem 19] stating that a necessary condition for supersolidity is the presence 
of either vacancies, or interstitials, or both. With this theorem, one just needs to make 
sure that there are groundstates where creating a vacancy and an interstitial cost finite 
energy. The latter is known to be the fact at least since Andreev and Lifshitz analysis 
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FIG. 11: Single-particle Green's function G(k = 0, r) computed by worm algorithm for hep 4 He 
at the melting density ra o =0.0287 A -3 and T = 0.2 K. Symbols refer to numerical data, solid 
lines are fits to the long-time exponential decay. The given numerical values are the interstitial 
(Ai = 22.8 ± 0.7 K) and the vacancy (Ay = 13.0 ± 0.5 K) activation energies, inferred from the 
slopes of G. The straight-line asymptotic behavior indicates that finite-temperature corrections 
are negligible. The inset shows the vacancy-interstitial gap -Egap — -|- Ay for different system 



sizes, proving that the results have reached their macroscopic limit. (Reproduced from Ref. 



2l|.) 



[20I ]: The energy for creating a vacancy/interstitial is positive, and arbitrarily close to the 
classical-crystal value, in the limit of strong interaction/large particle mass. 

The theorem of Ref. offers a reliable protocol of numerical proof that the groundstate 
of perfect 4 He hep crystal is an insulator. It is sufficient to demonstrate that (i) creating 
a vacancy/interstitial in the simulation box costs a finite energy, the value of which does 
not vanish with increasing the system size and decreasing the temperature, (ii) a state with 
finite concentration of vacancies/interstitials — that could potentially differ from the single- 
vacancy/interstitial situation due to collective effects — is unstable in the thermodynamical 
limit with respect to aggregation (i.e. a crystal at T = purges itself from the vacancies 
and interstitials). Both facts are demonstrated in Figs. [TI] and [12j (The aggregation of 
interstitials has not been studied since these cost much more energy than the vacancy, 
rendering the scenario of interstitial-induced supersolitity not realistic.) 

Having established that the perfect hep 4 He crystal is not supersolid, one has to explore 
disordered scenarios, when the superflow in a crystal is supported by defects. First-principles 
simulations performed by UMass-ETH-UAlberta-CUNY collaboration (briefly reviewed in 



Ref. 



181 ]) have revealed a number of disorder induced mechanisms of supersolidity in 4 He: 



18 



0.8 
0.6 
0.4 
0.2 





• • • • v a • • - 
• ••••• 



12 
r(A) 



15 



FIG. 12: The vacancy-vacancy correlation function v{r) as a function of the distance r between 
the vacancies shows that three vacancies easily cluster and form a tight bound state. The inset 
shows a typical snapshot of a layer of atomic positions averaged over the time interval [0, 0\ (filled 
black circles). It is seen that the three vacancies (triangles) have a tendency to cluster in layers. 



(Reproduced from Ref. 



2l|. 



superfluid dislocations, grain boundaries, ridges, and also a metastable amorphous super- 
solid, the so-called superglass. Here we confine ourselves with presenting the results for 
the superfluidity in the core of a screw dislocation — arguably the cleanest Luttinger liquid 
system in Nature. 

To visualize spatially inhomogeneous superfluidity in a worm algorithm simulation, one 
can employ two similar approaches. One approach is to plot the condensate density map. 
The other and more general approach (that also works for lower-dimensional systems with 
the genuine long-range order destroyed by fluctuations of phase) is to visualize the macro- 
scopic worldline loops responsible for non-zero winding numbers, and thus for the superfluid 
response. Identifying these loops in a given worldline configuration, projecting them from 
the (d + 1) dimensions onto a plane in the real space, and performing the average over a 
representative set of configurations, one obtains the winding-circle map of the superfluid 
region. The Luttinger liquid core of a screw dislocation visualized with this technique is 
shown in Fig. | 



19 




FIG. 13: Luttinger liquid in the core of the screw dislocation in solid He revealed by columnar 
winding-cycle map. Simulations correspond to the temperature 0.25 K and density 0.0287 A -3 . 
View is along the :r-axis in the basal plane — perpendicular to the core. Shown with large dots (in 
the lower half of the plot only) are the atomic positions in initial configuration. The unit of length 
is lA. (Reproduced from Ref. 22].) 



Problem of Deconfined Criticality. Flowgram Method 

The standard Ginzburg-Landau- Wilson (GLW) scenario of critical phenomena excludes 
generic continuous transitions between states which break different symmetries, thus imply- 
ing that the transition, if exists, is of the first order. An intriguing possibility of breaking 



GLW paradigm was proposed in Refs. 



23N25| for the so-called deconfined critical points 



(DCP) in two spatial dimensions. Nowadays the problem of DCP is one of the most ex- 
citing, and yet controversial topics in the theory of phase transitions. Remarkably, the 
field-theoretical model for deconfined criticality — to be referred below as the DCP action — 
is given by two identical complex-valued classical fields coupled to a gauge vector field, in 
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FIG. 14: Data collapse for the flows of the SU(2) symmetric DCP action. The line is a fit repre- 
senting the master curve. The horizontal axis is the scale reduced variable C(g)L. (Reproduced 
from Ref. Q-) 

three dimensions. Despite its simplicity and apparent closeness to the single-component 
counterpart (known to be in the inverted XY universality class of continuous phase transi- 
tions), the DCP action is not amenable to reliable analytic treatments because of its runaway 
renormalization flow to strong coupling at large scales. To establish the order of the phase 
transition in this model one has to resort to numerics. 

In Refs. [26, 127]], the order of the phase transition in the DCP action was studied by the 
worm algorithm (for the U(l)xU(l)- and SU(2)-symmetric actions, respectively). Within 
the given universality class, the optimal choice of microscopic model is the high-temperature 
expansion (cf. Sec. J of a (cubic-) lattice theory, resulting in the bond-current model with 
the following partition function (for the sake of definiteness, we present the answer for the 
SU(2)-symmetric case). 

Z = J2 2*itc Qbond expi-Hj), Hj = (g/2) ]T V %3 1$, (11) 

{J} i,j;a,b;fj,=l,2,3 

Here J^" is an integer non-negative bond current of the component a = 1,2, living on the 
bond (i,fJL). The bond subscript of a current is represented by the site, i, and direction 
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FIG. 15: Flowgrams for the short-range model Vij = g5ij. The lower horizontal line features 
the 0(4) universality scaling behavior; for g < g c ~ 0.95 all flows are attracted to this line. The 
upper horizontal line is the tricritical separatrix (marked as TP). Above it, flows diverge due to the 
first-order transition detected by the bi- modal distribution of energy. (Reproduced from Ref. j^.) 

/i = ±1, ±2, ±3 from this site. The direction-dependent subscript reflects the fact that for 
a given geometric bond there are two different bond currents of the same component a, and 
A> = ^ti ~Ji+n,-u * s their algebraic sum (here p, is a unit translation vector in the direction 
/i). The bond currents are subject to the conservation constraint on each site: YIuH^j, = 0- 
The parameter t controls the strength of the lattice gradient term for the complex fields in 
the DCP action, and g is the coupling constant for the interaction between the complex and 
gauge fields. The integration over the gauge field results in the long-range interaction, Vy, 
between the currents. The Fourier transform of is given by V" q = 1/ ^^ =1 2 3 sin 2 (g M /2) 
and implies an asymptotic behavior Vy ~ 1/Vjj at large distances. It is this Coulomb 
asymptotic tail of the current-current interaction that leads to a qualitative difference of the 
DCP action from its short-range counterparts. 

With its closed-loop structure enforced by the current conservation constraint and the 
positive-definitness of the weighting factors, the model (TTTT) is in the domain of applicability 
of the worm algorithm, each of the two components being updated by its individual pair of 
worms. The results of simulations of both U(l)xU(l)- and SU(2)-symmetric actions lead 
to an unfortunate for the DCP theory conclusion that the phase transition is of the first 



22 



order. The definitive conclusion is based on the flowgram method [26] of finite-size analysis. 
The key idea is to demonstrate that the universal large-scale behavior at g — > is identical 
to that at some finite coupling g = g coll where the nature of the transition can be easily 
revealed. The procedure is: 

(i) Introduce a definition of the critical point for a finite-size system of linear size L 
consistent with the thermodynamic limit and insensitive to the order of the transition. 
Specifically, for any given g and L the critical value of t was defined by the requirement 
that the ratio of statistical weights of configurations with and without windings be 
equal to a fixed constant. 

(ii) At the critical point, calculate a quantity R(L, g) that is supposed to be scale-invariant 
for a continuous phase transition in question, vanish in one of the phases and diverge 
in the other. Specifically, once can take R(L, g) to be the variance of the winding 
number in the counter-flow channel: R(L,g) = (W 2 ) = — W^) 2 ), where 
W a ^ is the winding number of the component a in the direction /z. (See Ref. for 
the motivation of this choice.) 

(iii) Perform a data collapse for flowgrams of R(L,g), by rescaling the linear system size, 
L — > C(g)L, where C(g) is a smooth and monotonically increasing function of the 
coupling constant g. [In the present case, it is a priori known that C(g — > 0) oc g.} 

A collapse of the rescaled flows within an interval g G [0, g co \i] implies that the type of the 
transition within the interval remains the same, and thus can be inferred by dealing with 
the g = g co n point only. Since the g — > limit implies large spatial scales, and, therefore, 
model-independent runaway renormalization flow pattern, the conclusions are universal. 

As is seen in Fig. [TH the flows for the DCP action collapse perfectly in the region 
0.125 < g < g co n = 1.65. The rescaling function C(g) exhibits the expected linear behavior 
C(g) oc g at small g. 

In accordance with the above-outlined logic of flowgram method, the flow collapse within 
the interval g G [0, f? co ii] proves that the order of the transition within this interval does 
not change. The mere fact of the data collapse on a master curve with a finite slope 
is not sufficient to conclude that the transition is of the first order. What appears to 
be a (characteristic of the first-order transition) diverging behavior in Fig. [TH might be 
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just a reconstruction of the flow from the 0(4)-universality (at g = 0) to a novel DCP- 
universality at strong coupling. To complete the proof, one has to determine the nature of 
the transition for g = g coll . In this parameter range, the standard technique of detecting 
discontinuous transitions by the bi- modal energy distribution becomes feasible. At g — 1.65 
it becomes possible to clearly see a bi-modal energy distribution that gets more and more 
pronounced with increasing the systems size. This leads to a conclusion that the whole 
phase transition line for small g features a generic weak first-order transition. Driven by 
long-range interactions, this behavior develops on length scales oc 1/g — > oo for small g and 
thus is universal. 

Finally, it is very instructive to contrast the flowgram for the DCP action (TiTj) with the 
flowgram for the short-range counterpart of (fTTj) , where = g5ij. The flows for the short- 
range model are presented in Fig. [151 As opposed to the DCP action case, it is impossible 
to collapse the data on a single master curve by re-scaling L. Now the flows clearly reveal 
a tri-critical point separating the second-order part of the phase transition line from the 
first-order part. 



CONCLUSIONS AND OUTLOOK 



Worm Algorithm is a technique for performing efficient updates of configurations that 
have the form of closed paths/loops. The closed-loop structure imposes topological con- 
straints and thus creates ergodicity problems for local updates in large system sizes. Worm 
Algorithm works in an extended configuration space containing all the original configura- 
tions as well as configurations with open loops. All updates are local and are performed 
exclusively at the open loop ends, referred to as worms. In most cases, the open-loop configu- 
rations themselves are of prime physical interest being associated with correlation functions, 
such as the single-particle Green function. It was demonstrated for a variety of systems and 
universality classes that WA eliminates problems with ergodicity and critical slowing down, 
successfully competing even with model-specific cluster algorithms 28j. At the same time, 
it is a flexible approach with a broad range of applications in many-particle bosonic and 
spin systems. It readily produces loops winding around the system, allows efficient simu- 
lations of off-diagonal correlations, grand canonical ensembles, disordered systems, etc. At 
the moment WA has no competitors among unbiased first-principles approaches for bosons 
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with strong interactions between the particles, as well as with strong external — both regular 
and disordered — potential. It is easy to predict that in the nearest future WA will remain 
the method of choice for detailed studies of non-trivial strongly correlated bosonic systems 
(multicomponent, disordered, with long-range interactions, solid and supersolid 4 He, etc.) 
Recently, WA has proved indispensable for guiding experimental efforts in creating optical 
lattice emulators; it will continue playing this important role. 

In its most general form, the idea of WA is to work in an enlarged configuration space 
which includes configurations violating constraints characteristic of the physical configura- 
tions. Nowadays worm-type updates and/or worm-type estimators for the Green function 
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32]. 



are an integral part of many other state-of-the-art lattice Monte Carlo algorithms 
It is only upon the implementation of the worm-type updates it became possible to over- 
come the critical slowing down in the Stochastic Series Expansion scheme 311 ] . Important 



improvements have been made to achieve maximal efficiency of the worm-type update s by 
suppressing bouncing: The so-called directed (guided) loop and geometrical WA [3l|, l33| . 
appear to be an optimal combination in terms of universality and performance. Successful 
applications of WA in high-energy physics 34M36I] . with exciting most recent developments 
render the approach interdisciplinary. Whenever a new representation for Quantum 



appncc 
37H39 



Monte Carlo appears, the generalized WA idea may prove useful for developing an efficient 
updating strategy. 
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